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The scattering of electromagnetic wave by a periodic array of nanowires is cal- 



culated by the boundary element method. The method is extended to the infinite 



grating near the interface between two dielectrics. A special Green function is de- 
rived that allows to study the evanescent wave. The Rayleigh — Wood's anomalies 
are found in the period-to-wavelength dependence of the average Pointing vector in 
the wave zone. For thin wires the calculations are shown to agree with the two- 



dimensional coupled dipole approximation. 
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I. INTRODUCTION 

The travelling wave falling onto the boundary between two dielectrics by the angle of total 
internal reflection generates an evanescent wave in the second medium. The illumination 
by evanescent wave is widely used in applications. In near-field scanning optical microscopy 
it allows one to detect the desired signal with no background Total internal reflection 
fluorescence, exploiting decay of evanescent wave, achieves nanometer resolution along the 
optical path through a thin object |2J. Covering the fiber surface by nanowires improves the 
sensitivity of fiber Bragg grating refractometry 

While the scattering problem for evanescent wave is of importance, it has, however, no 
analytical solution. The necessity of treating the evanescent wave along with its source 
exaggerates the calculation. For travelling wave the analytical solution exists for the homo- 



geneous medium and one circular cylinder 



Q. 



The numerical methods allows one to solve 



the problem for one elliptic cylinder [5J. The numerical solution for two different dielectric 



half spaces are found recently for one [6j and two [7] cylinders. The scattering by few cylin- 
ders has very small cross-section and then is difficult to observe. To enhance the effect one 
needs to consider a number of cylinders, but growing the number of scatterers results in 
additional computational difficulty To remedy this contradiction between the observation 
needs and computation difficulties, we reduce a two-dimensional infinite periodic array to 
only one unit cell with the help of periodicity. The aim of the present paper is to calculate 
the scattering of evanescent wave by the periodic array of parallel cylinders of arbitrary 
diameter. 

The periodic structures of scatterers with subwavelength dimension of each enhance the 



signal compared to one scatterer 8KH|. Two-dimensional structures are important, for 



instance, for photonic crystals lfj or metamaterials [l2|. The diffraction grating can be 
characterized by its Rayleigh — Wood's anomalies. We consider these resonances and com- 
pare their structure with the case of travelling wave and with widely used coupled dipole 
approximation. 

For this purpose we propose a modification of the boundary integral equations, Section 



HT1 which are the basis of the Boundary Element Method (BEM) 13H15J. We derive a special 
Green function that takes into account the evanescent wave together with its source, the 
dielectric half-space, in Subsection III Al The method is extended to the periodic array in 
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Subsection III Bl In SectionQTTlwe present the results of calculation. Subsection IIIIAI contains 
the study of average energy flux in the far-field domain as a function of period-to- wavelength 
ratio, including splitting of Rayleigh — Wood's anomalies. The Subsection IIIIBI is devoted 
to the calculation within dipole approximation, in comparison to the more rigorous BEM, 
for the travelling exciting wave. The coupled dipoles equations are presented for both finite 
and infinite periodic grating. 



II. METHOD 



A. Green Function 



Let us consider two media, 1 and 2, separated by a plane boundary at y — 0. We look 
for the specific Green function G(x, y; x', y 1 ) satisfying inhomogeneous equation 



AG+ k 2 G = 8(x - x')5(y - y') 



(1) 



in these media. Function G depends on the difference x — x' only due to translational 
symmetry. The boundary condition at y — for p-wave is 



[G(x, y; x', y')] 



After the Fourier transformation 



y=o 



ldG{x,y;x',y') 
e 



dy 



0. 



(2) 



1 f 00 

G(x,y;0,y') = — J G g (y,y')e^dq 



(3) 



Eq. (IT]) is reduced to an ordinary differential equation having exponential solutions. Using 
conditions ([2]) we obtain the function. The expression for y' > in g-domain is 



G n 



[1 + r(q)) e 



y<o, 



-w\y-y'\ _)_ r (g) e -M2(y+?/') ) y > 0. 



(4) 



Here 



r(q) 



£lfl 2 + 

is the Fresnel reflection coefficient of p-wave at normal incidence 



16|, = q 2 



(5) 
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Carrying out Fourier transformation ([3]) of function t^j at y > we have two terms 

oo— iO 

C -M2 [if-V I +i«(a:-a:') ^ 



— oo+iO 
oo— iO 



-oo+iO 



47T/i 2 

^2(y+y')+ig(x-x') r ( ( l) d 1 
47T//2 



(6) 



(7) 



where G = Gi + G 2 . The sign of the square root is given by the rule y / g 2 — k\ — > 
-i^kl-q 2 ^ 2 < k\. 

The first term ([6]) can be calculated analytically and reduces to the known Green function 



of Helmholtz equation in the homogeneous space 



GAr-r 1 



(hp. 



4i 



where denotes the Hankel function of the first kind [1 



p ± = ^(x-x'Y + (y±y'Y 



(8) 



(9) 



are the distances from the source and its image. The second term G 2 gives the effect of 
the reflected image source. The amplitude of the source at each q is equal to the reflection 
coefficient r(q). Thus along with the point source at [x',y') we have to consider mirror- 
image source r(q) at (x', — y'). The total field at the upper half-plane is the sum of the fields 
generated by the two sources at each q. This approach intrinsically takes into account the 



multiple scattering. The Green function of this type was studied 



spherical acoustic, see 



19] , or cylindrical electromagnetic waves 



'or homogeneous waves: 



20] . A general approach 



has been built, based on integral equation, to study the field scattered by an object near 



two-media interface, for excitation by homogeneous 



2l| and evanescent 



flflli 



ight wave. 



Analogous calculations give the Green function for the source inside the lower dielectric. 



G a 



2/n 



e -m\v-y'\ _ r (q)e fll{y+y '\ y < 0, 
[1 - r(q)]e~L l2y+ ^y\ y > 0. 



(10) 



Comparing function f flQ|) with @ we see that the expression can be obtained by the sym- 
metry transformation, namely, changing medium indices 1 <H- 2 and the sign of the ordinate 
y —7- —y,y f —> —y'. The case of s-wave the boundary conditions (T2]) are replaced by 

~dG(x,y;x',y') 



[G(x, y; x , y% =0 



dy 



0. 



(11) 
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FIG. 1. A cross-section of the periodic array of cylinders: a is the radius, h is the gap, I is the 
period, 9\ is the light incidence angle. 

The result can be also obtained from (jlj) replacing the Fresnel reflection coefficient (jSJ) for 
normal incidence by its expression for s-wave 

fi 2 - Hi 



r(q) 



(12) 



Hi + H2 

The compound Green function is valid in both half-spaces and makes it possible to solve 
the scattering problem near the interface between two dielectrics. This function is yet for 
few cylinders but worthless for the periodic grid. In the next section we explain how to 
extend the computation for the infinite grid with the help of Floquet theorem and derive 
the corresponding effective Green function. 



B. Boundary Element Method 

Let us consider a scattering of a plane wave by the periodic grid of identical parallel 
cylinders with dielectric permittivity £3, see Fig. [U The grid is located in dielectric half- 
space 2 (with dielectric permittivity £2) parallel to its boundary with half-space 1 (si). A 
travelling wave with frequency u falls from half-space 1 by the angle of incidence 6\ . Denote 
&i,2 — &o^/ £ i,2 the wavenumber in each dielectric, where k = u/c, c is the speed of light. As 
known from Fresnel formulas, the transmitted wave in media 2 becomes evanescent when 
angle B\ exceeds the angle of total internal reflection. Choose axis z along the axes of 
cylinders, axis x along the boundary, and axis y perpendicular to the boundary, as shown 
in Fig. [TJ 

We concentrate on the p-wave (TM case), when the magnetic field has only ^-component, 
which consists of the falling and the scattered fields: H z = H + H'. With the help of integral 
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Green formula we write the coupled boundary integral equations for scattered field u = H' 
and its normal derivative / = dH' /On at a contour point P = (x, y) G Y: 



u(P) 



+ 



dG 



out 



. dn Q 



u{Q) - G out f(Q) 



dQ = 0. 



(13) 



Here riq is the internal normal in point Q = (x',y') G T, G out (P,Q) is the external Green 
function. The integration point Q runs through compound contour F = U m T m , where r m is 
the contour of m-th circle (cross section of cylinder): m — 1, 2, . . . , N; N is the number of 
cylinders. The integral over an unlinked domain is a sum of integrals over partial contours 
T m . This equation corresponds to the outer limit, i.e. the limit of Green formula when the 
observation point tends to contour T from the external side. For the outer limit we use the 
special Green function described in Sec. Ill AI It explicitly takes into account the boundary 
conditions at theplane interface between the dielectrics 1 and 2, and has a form of divergent 
wave at infinity |6|, |7|. 

The second integral equation arises as a result of the passage to the inner limit, i.e. when 
point P tends to the contour from internal area: 



u(P) 



~ dG m 
_dn Q 

Uq _ 
2 



u(Q) - G m f(Q) 

dGi n 

-Uq 



dn 



dQ = B , 
Ginfo dQ. 



(14) 
(15) 



Here G in (P, Q) is the internal Green function, uq = Ho, fo = dHo/dn. For internal domain 
we use the cylindrical wave 

1 



G- m (x,y;x',y') = ^-U^ifap 



(16) 



For iV cylinders there are N internal and iV external equations. Solving the set of 2N 
equations we find functions u(P),f(P),P G T. If we know functions u, f at the contour, 
then scattered field in arbitrary external point R = (x,y) ^ T can be found from the Green 
theorem 

dG out 



u(R) 



-u(Q) - G ou J(Q) 



dQ, 



(17) 



where G out = G out (R, Q) is the external Green function. 

Eqs. ( 1T31) . ( IT4l) refer to the case of a finite number iV of the cylinders. Passing to 
the limit iV — > oo we consider points at the compound contour with equal y m = y and 
Xm — XQ + ml;m = 0, ±1, ±2, ... A zeroth cylinder m = is chosen arbitrary for the infinite 
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grating. The following property (also known as Floquet theorem) is a consequence of the 
periodicity: 

u(x m , y m ) = u(x + md, y ) = u(x , y )e lk±ml , (18) 

where k± = k\ sin 9\ is the tangential component of the wavevector conserved at the bound- 
ary. Taking into account the same property of function f(Q) and reversing the order of 
summation and integration in Eq. ([TBI and ffT7|) 22] we obtain the sum: 



E 



m=— oo ■- 



dG on t 
dn 



u{Q) - G ont f(Q) 



Ak i ml 



(19) 



The special Green function ([3]), (jlj) in the upper half space 2 can be written as a contour 
integral 



Gr 



e -Hw-VI + rke 



c 



(20) 



where y, y' > 0, A = x — x', G out = G out (x, y; x', y'), /i = \/k 2 — k%, C = (— oo + zO, +oo — iO). 
Consider the latter sum in Eq. (|T9|) 

+oo 



G cS = J2 G out (x,y;x' m ,y' m y k - ml . (21) 

m=— oo 



From ( )2T1) and integral representation (l20l) . alternating the order of summation and integra- 
tion, we obtain the periodic delta function 

+00 



E 



m=— 00 



-ikmd+ik ±md 



21T 

—5 {k — k±_ — nx) 

v 



(22) 



where n = 0, ±1, ±2, . . . , x = 2n/l is the reciprocal lattice vector. Then integrating the 
delta-function over k we get the expression for G c fr 



G 



+°° p i(fe±+>m)A 



eff 



E 



2//x„ 



g-Mn|2/-2/ol _|_ rn g-Mn(2/+2/o) 



(23) 



where A = x — a;' , /i n = a/ (k± + nx) 2 — fcf, r n = r fc±+n ^. The former term in Eq. f|T9l is 
reduced to the normal derivative of ( 1231) analogously, changing the order of normal derivative 
and summation. 

Expression ( |23l) determines the effective Green function, taking into account the infinite 
grid as a whole. If we replace G out by G c g in external integral equation (fTB"]) and Green 
theorem ffTTl) . then the path of integration will be r . It is the zeroth circle, the connected 
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contour, instead of infinite unlinked contour T. Only one internal equation should be solved, 
for m = 0. After finding functions u(P), f(P) we use Eq. (fTT|) for calculation of the field in 
the wave zone. 



III. RESULTS 
A. Rayleigh — Wood's anomalies 

We study the transparency of the grating as a function of distance I between the neighbor 
cylinders. The transparency is defined as a ratio of y-projection of average Pointing vector to 
the energy flux density of the falling wave. The Pointing vector is averaged over the period 
and calculated at sufficiently large distance from the boundary y ~ 10A, where A = 1.5 /im 
is the fixed wavelength. 

Only terms corresponding to travelling waves with (k± + nx) 2 < k\ (and then Re/i ra = 0) 
have to be kept for calculations in the wave zone. Thus, the sum (123]) becomes finite and 
relatively easy for numerical calculation compared to integral representation (|20|) . Fig. [2] 
shows the transparency as a function of the spacing between the cylinders, I, (rated to the 
wavelength) for the travelling wave (ei = e 2 = = 45°). The grating is placed close to 
the surface: h = 0.01 //m. 

We see the narrow dips, so-called Rayleigh — Wood's anomalies (RWA). These resonances 
take place when the grazing mode appears. The condition of resonance is synchronization 
of the phases between the waves scattered by neighbor cylinders. The grazing modes may 
have both positive and negative order n = ±1, ±2, . . . The positive orders n > generate 
the surface wave going in negative direction and vice versa. For positive /negative mode the 
synchronism condition is 

k 1 lsm6 1 ±k 2 l = 2n\n\. (24) 

Here /ci/sin^ is the difference of optical paths in half-space 1, k 2 l is the difference in half- 
space 2. The expression (1241 is valid only for small diameters 2a <C A. 

Let us express the wavenumbers in terms of the wavelength ki >2 = 27r v /£i i2 /A. Eq. ( )24l) 
can be rewritten as 

I ±\n\ 
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FIG. 2. Raleigh — Wood's anomalies in the transparency for travelling wave {e\ = 1, the cylinders' 
radii a = 0.1 /iin). The first dip with higher resolution is shown in the inset. 

In vacuum for n = 1,2 and 9\ = 45° it gives l/X = 0.5858, 1.1716, . . . Near the resonance 
a rearrangement of the wave pattern occurs due to arising the next order diffraction and 
transparency decreases sharply. The inset demonstrates the splitting of the first dip. Its 
physical nature is a difference between induced dipole moments along x and y directions. 
This difference is verified in the section IIIIBI within the coupled dipole approximation. 

The case of evanescent wave is shown in Fig. [31 The evanescent wave in the upper 
half space 2 arises from the travelling wave falling onto the boundary by the angle of total 
internal reflection (TIR). The refractive indices are y/ei = 1.5, y/e2 = 1 (describing a glass 
in the air), the incident angle Q\ = 45°, as above. The limit value of transparency at I ^> A 
is zero instead of unity. The positions of resonances have shifted, since in dielectric Eq. ( )25l) 
gives l/\ = 0.4853, 0.9706, . . . 

Comparing Fig. [2] and [3] we see an essential difference. There are four main distinctions 
of Fig. |3j the RWA broadening, the vanishing splitting, the weak step ("jag") at l/X w 0.78, 
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FIG. 3. The transparency as a function of period-to-wavelength ratio for evanescent wave (ei 
2.25). An arrow indicates the weak step. 



and the zero asymptotic value at l/X — > oo. There are two differences in conditions between 
Fig. [2]and[3j the dielectric lower half-space with ^fe{ = 1.5 and the TIR incident angle. To 
clarify the contributions of both differences let us repeat the calculation for two very close 
angles, less and more than the TIR angle. We fix e\ = 2.25, 62 = 1 and choose two values 
in incidence angle 9\, both very close to the angle of total internal reflection. For our case 
the TIR angle is equal to #0 — arcsin (l/v/el) = 0.729728. Our values are 07/ « 0.7297 and 
61 ~ 0.7296. First is less than TIR angle, the second is greater. The difference is very small, 
\6f-8 \ « 10- 4 . At the same time the angles corresponds the travelling or evanescent wave 
in medium 2 at almost equal conditions. The transparency curves are shown in Fig. HI The 
upper curve, the transparency for travelling wave, differs noticeably from the lower one. 

When the angle 6\ is close to the TIR angle 6q, only RWA of positive order exist at 
l/X = n/2,n = 1,2,... The travelling wave propagates nearly along the surface at re- 
direction and has mainly ^-component of the electric field. Then the x-component of the 
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FIG. 4. The transparency as a function of period-to-wavelength ratio for a narrow gap h = 0.01 //m: 
the incidence angle is 9f — evanescent wave (curve 1) or 9± — traveling wave (curve 2). An arrow 
indicates the weak step. 



dipole moment is very small and the splitting vanishes. The same situation takes place for 
angles 9\ > 9q, i.e. for evanescent wave, when x-component of the dipole moment is exactly 
zero. Unlike, at a smaller angle 9\ the x-component of dipole moment arises along with the 
RWA splitting. Also, the broadening takes place for both the curves. 

The picture considerably varies also with changing the gap h. The computation in Fig. H] 
correspond to h = 10 nm. Fig. shows the same curves for the grating placed far from the 
boundary, at h — 300 nm, when the role of multiple reflection becomes relatively weak. In 
comparison with Fig. H] the gap increasing results in higher peak values of the transparency 
coefficient. The step at l/X ~ 0.78 vanishes for wide gap, and then it can be interpreted as 
an effect of the image source. 

Thus the broadening and one component RWA are the effects of dielectric. The jag is 
the result of image grating reflected in the dielectric. The zero limiting value at l/X — > oo 
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B. Discrete Dipole Approximation 



In this section we consider the grid of nanowires within dipole approximation. It means 
that we will describe the response of the wire on the exciting electric field in terms of the 
wire's dipole polarizability, a. As before, we consider the case when the electric field lies in 
xy-plane, the acquired dipole moments are 2-dimensional. Also, the dipole-dipole interaction 
between the nanowires should be taken into account. Such an approach dealing with coupled 
dipoles is widely used in nanophotonics, particularly, in optics of nanoparticles 

Qfl- The 

concept, basically, consists in finding the values of the dipole moments, acquired by the 
particles, when each is excited by both the incident (given) field and the sum of the dipole 
fields from all remaining particles. The dipoles obey the so-called Coupled Dipoles Equation 
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(CDE): 

ar 1 d i = 5^G y d i + E w , (26) 

where Eo« = Eo(rj) — the electric field of the incident electromagnetic wave at the point 
of z-th dipole's location; «j is the polarizability of i-th particle (the nanowire, in our case); 
summation is made over the whole ensemble of the wires. The diadic Gij is defined through 
the electric field produced by j-th dipole at point iv 

E a ( r i) = Gij >a pdj t p, (27) 

where Greek indices denote Cartesian components. Thus, Gijdj is the contribution from 
j-th dipole to the field, acting on i-th one. 

After the system of CDE has been solved, any quantity of interest can be found using 
the values of dj. Here we calculate the electric field near the nanowares array, to compare 
with the results produced by BEM method. 

In Figj6] the squared modulus of the electric field, calculated at the center of the chain of 
500 gold nanowires, is plotted for various spacing between them, I, and their radii, a. The 
exciting field was the plane wave falling by the angle Q\ = 7r/4 from the left-bottom side (Fig. 
[H the polarization is in xy plane; the wavelength is A = 1512 nm. As we see, the results of 
two compared methods perfectly match for 30-nm and 50-nm cylinders, while for the case 
of 100-nm radius there is a noticeable difference. The discrepancy arise when the diameter, 
2a, becomes comparable with \/2it, destroying the dipole approximation. Quite remarkable 
is that such a good agreement between two methods remains up to the gap between the 
cylinders being as small as their radius. 

The computational difficulty of CDE solving rapidly grows with the number of interacting 
dipoles. However, there is a special case - - an infinite chain of dipoles, for which the 
calculations become drastically more simple. It allows to consider the infinite grating of 
nanowires analytically, at least for excitation by plane travelling wave. 

The main idea is that if the grating of nanowires and the aperture of the incoming wave 
are both infinite, the wires become fully equivalent. So, the difference in the corresponding 
dipole moments may originate only from the different values of the incident field at each 
dipole. In other words, the contribution to the total (acting) field, which is due to dipole- 
dipole interaction, J2j=£iGijdj, must have the same amplitude at each i-th dipole, and must 
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FIG. 6. Comparison of CDE solution for the chain of 500 2D dipoles (solid lines) with that obtained 
by BEM for the infinite array of nanocylinders (circles and squares). Plotted is the value of squared 
electric field in the middle of the gap between the cylinders (for CDE, between 250-th and 251-th 
ones). Here a = 30 (curve 1), 50 (2), and 100 nm (3). The field is normalized by the amplitude of 
the falling wave, Eq. 

be always in phase with the incoming field Eo(rj). Applying these two conditions to f )26|) 
we can get an analytical solution of CDE. For this purpose, we use the explicit expression 
for the electric field of an oscillating 2-dimensional dipole moment d oc exp (—iut) 



Efr) 



ink 



r 2 d (krH ( o\kr) - H^\kr) 



r (d ■ r) ( krU^ {kr) - 2U\ l) (kr 



(28) 



which is derived using the approach, very similar to that described in 27| for 3-dimensional 
case. Also, we use the asymptotic expression for both Hankel functions, that gives an 



15 



inaccuracy less than 15% for kr > 0.5. After summation we get: 



d x = E Q a cos 6i x 



(29) 




i 



d y = E a sin 9\ x 



(30) 



(1 + i)y/na 



i 



1 



[8kH 2 F 1/2 + 7iklF 3/2 - 3F 5/2 ] 



where 



F s = Li s (e ikl{1 ~ sin9l) ) + Li s (e ikl{1+Sin9l) ), 



lA s (z) is the polylogarithm [17] 



n=l 



The dipole polarizability, a, is the same for all the nanowires; k = u/c. 

The resulting dipole moment d x is plotted in FigJJJ The calculations were done for gold 
nanowires with the radius 100 nm as a dependence of the nanowires spacing, at constant 
wavelength of incoming radiation A = 1512 nm. The Rayleigh- Wood's anomalies are clearly 
seen in Fig. [7J as peaked extrema. 

For comparison, we performed also a direct numerical solution of CDE (|26|) for 500 two- 
dimensional dipoles, formed by the gold nanowires. The results for d x component (taken for 
the "middle" dipole, i = 250) are in excellent quantitative agreement with the analytical 
given above, as also shown in Fig. [71 

The results of the same calculations for d y are shown in Fig. [8] Comparing to Fig. [7J we 
see that positions of peaks do not coincide with the case of d x . It is the shift that produces the 
splitting shown in Fig. [2j As for the amplitudes difference between analytical and numerical 
solutions, we believe it to be a consequence of high derivatives near the resonances, that 
drops down the numerical accuracy. 



An evanescent wave has to be considered together with its source. The special com- 
pound Green function is proposed taking into account the falling waves in the lower half 
space. Boundary integral equations are extended for periodic infinite structure. The effec- 
tive Green function is found allowing to treat an infinite periodic array of nanowires. The 



IV. CONCLUSIONS 
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FIG. 7. The absolute value of x-components of the dipole moment d x vs the relative spacing 
between the wires, l/X: the analytical formula for infinite number of wires (solid line); numerical 
calculation for array of 500 elements (dots). 



numerical calculation for circular cylinders reveals Rayleigh - - Wood's anomalies in the 
grating transparency. The dips turn into peaks and splitting of resonances vanishes when 
total internal reflection occurs. For the case of evanescent wave nulling the asymptotic value 
of transparency is shown. The calculations are compared with coupled dipole approximation 
using the analytical solution, derived for an infinite number of coupled dipoles, and the nu- 
merical (direct) calculation for the finite grating with N = 500 wires. In the latter case the 
oscillations are found in the near-field, which are the result of a finite number of cylinders. 
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